# Dickstein, Ho, and Mark (2023)
# This script analyzes how the mandated insurance counterfactual
# affects consumer surplus for subsets of the sample.

# ~ # ~ # ~ # ~ # 
# Preliminaries #
# ~ # ~ # ~ # ~ # 
setwd("../../../../library")
source("PreliminariesCode.R")

# What alpha-psi cutoff to use: 
alphampsicutoff <- .05

# What ACG cutoff to use: 
acgHL_boundary <- 1.76

# What is nu? 
nu <- .65

# what is the base markups to use?
mrkup_vec <- c(0, 25)

# Which spec do you want to use?
spectouse <- "lowrisk_simplemh_censor0.2"

# Where is the counterfactual folder?
counterfactual_folder <- paste0(project_folder, "/analysis/counterfactuals")

# Loading counterfactual functions 
source(paste0(counterfactual_folder, "/library/DA03aRunCounterfactualFunctions.R"))

# Create directory to save in
dir.create(paste0(project_folder, "/analysis/tablesandfigures/release/misc"), recursive = T)

# ~ # ~ #
# Data  #
# ~ # ~ #

# Load the household data
finalsubsdata <- fread(paste0(project_folder, "/data/all_markets_subscribers.csv"))

# Split the data up by type
finalsubsdata$fpl_quantiles <- cut(
  finalsubsdata$best_guess_incomeoverFPL, 
  quantile(finalsubsdata$best_guess_incomeoverFPL, seq(0, 1, .25)))
fpl_quantile_vec <- unique(finalsubsdata$fpl_quantiles)
fpl_quantile_vec <- fpl_quantile_vec[which(!is.na(fpl_quantile_vec))]
fpl_quantile_vec <- fpl_quantile_vec[order(fpl_quantile_vec)]
finalsubsdata$risk_quantiles <- cut(
  finalsubsdata$sum_concurrent_risk, 
  quantile(finalsubsdata$sum_concurrent_risk, seq(0, 1, .25)))
risk_quantile_vec <- unique(finalsubsdata$risk_quantiles)
risk_quantile_vec <- risk_quantile_vec[which(!is.na(risk_quantile_vec))]
risk_quantile_vec <- risk_quantile_vec[order(risk_quantile_vec)]
finalsubsdata[, ftype := fcase(
  withkids == T & age > 50, 1, 
  withkids == F & age > 50, 2, 
  withkids == T & age <= 50, 3, 
  withkids == F & age <= 50, 4
), ]
ftype_vec <- unique(finalsubsdata$ftype)
ftype_vec <- ftype_vec[order(ftype_vec)]


for(mrkup in mrkup_vec){
  # Load the base data
  load(paste0(counterfactual_folder, "/specs/", spectouse, "/output/OnlyIndividualMarketCounterfactual_mrkup", mrkup, ".RData"))
  load(paste0(counterfactual_folder, "/specs/", spectouse, "/output/CombinedMarketCounterfactual1_mrkup", mrkup, ".RData"))
  load(paste0(counterfactual_folder, "/specs/", spectouse, "/output/CombinedMarketCounterfactual7_mrkup", mrkup, ".RData"))

  # Load the SG surplus data
  sg_surplus <- fread(paste0(counterfactual_folder, "/specs/", spectouse, "/data/SGSurplusEstimationFile.csv"))
  ## Dropping outside options and keeping 2016
  sg_surplus <- sg_surplus[!(x_av == 0)]
  sg_surplus <- sg_surplus[year == 2016]
  sg_surplus <- sg_surplus[!(fixedeffect_beta0 == 0 & x_av != 0)]

  ## Collapse to subscriber ID: 
  sg_surplus <- sg_surplus[x_av > 0, 
    .(logsum = log(sum(exp(V))), alphampsi = mode1(alpha - psi)), 
    by = c("subscriberid", "year")
    ]
  sg_surplus$cs <- sg_surplus$logsum / sg_surplus$alphampsi

  # Calculating Results for different population subsets:
  ## Calculating results for overall
  overall_resbase_list <- list()
  overall_res5_list <- list()
  print(paste("Starting overall"))
  overall_resbase_list[[1]] <- generateOutcomeColumns(
    CounterFactualResults = list(
    ExplodedData = OnlyIndividualMarketCounterfactual$ExplodedData,
      PlanData = OnlyIndividualMarketCounterfactual$PlanData),
    SGpremsubs = T, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio")
  overall_resbase_list[[1]][[2]] <- generateSGOutcomeColumnsPre(
    CounterFactualResults = list(
      ExplodedData = CombinedMarketCounterfactual1$ExplodedData,
      PlanData = CombinedMarketCounterfactual1$PlanData),
    SGpremsubs = T, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio")
  overall_resbase_list[[1]][[2]]$welfare[1] <- mean(sg_surplus[alphampsi > 0 ]$cs)
  overall_resbase_list[[1]][[2]]$welfare[2] <- mean(sg_surplus[alphampsi > alphampsicutoff]$cs)
  overall_res5_list[[1]] <- generateOutcomeColumns(
    CounterFactualResults = list(
      ExplodedData = CombinedMarketCounterfactual7$ExplodedData,
      PlanData = CombinedMarketCounterfactual7$PlanData),
    SGpremsubs = F, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Fixed")
  
  ## Calculating results for income bins
  fpl_resbase_list <- list()
  fpl_res5_list <- list()
  for(i in 1:4){
    print(paste("Starting quantile vec", i))
    fpl_resbase_list[[i]] <- generateOutcomeColumns(
      CounterFactualResults = list(
        ExplodedData = OnlyIndividualMarketCounterfactual$ExplodedData[orig_subs_id %in% finalsubsdata[fpl_quantiles == fpl_quantile_vec[i]]$subscriberid],
        PlanData = OnlyIndividualMarketCounterfactual$PlanData),
      SGpremsubs = T, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Ratio")
    fpl_resbase_list[[i]][[2]] <- generateSGOutcomeColumnsPre(
      CounterFactualResults = list(
        ExplodedData = CombinedMarketCounterfactual1$ExplodedData[orig_subs_id %in% finalsubsdata[fpl_quantiles == fpl_quantile_vec[i]]$subscriberid],
        PlanData = CombinedMarketCounterfactual1$PlanData),
      SGpremsubs = T, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Ratio")
    fpl_resbase_list[[i]][[2]]$welfare[1] <- mean(sg_surplus[alphampsi > 0 & subscriberid %in% finalsubsdata[fpl_quantiles == fpl_quantile_vec[i]]$subscriberid]$cs)
    fpl_resbase_list[[i]][[2]]$welfare[2] <- mean(sg_surplus[alphampsi > alphampsicutoff& subscriberid %in% finalsubsdata[fpl_quantiles == fpl_quantile_vec[i]]$subscriberid]$cs)
    fpl_res5_list[[i]] <- generateOutcomeColumns(
      CounterFactualResults = list(
        ExplodedData = CombinedMarketCounterfactual7$ExplodedData[orig_subs_id %in% finalsubsdata[fpl_quantiles == fpl_quantile_vec[i]]$subscriberid],
        PlanData = CombinedMarketCounterfactual7$PlanData),
      SGpremsubs = F, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Fixed")
  }
  
  ## Calculating results for risk bins
  risk_resbase_list <- list()
  risk_res5_list <- list()
  for(i in 1:4){
    print(paste("Starting quantile vec", i))
    risk_resbase_list[[i]] <- generateOutcomeColumns(
      CounterFactualResults = list(
        ExplodedData = OnlyIndividualMarketCounterfactual$ExplodedData[orig_subs_id %in% finalsubsdata[risk_quantiles == risk_quantile_vec[i]]$subscriberid],
        PlanData = OnlyIndividualMarketCounterfactual$PlanData),
      SGpremsubs = T, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Ratio")
    risk_resbase_list[[i]][[2]] <- generateSGOutcomeColumnsPre(
      CounterFactualResults = list(
        ExplodedData = CombinedMarketCounterfactual1$ExplodedData[orig_subs_id %in% finalsubsdata[risk_quantiles == risk_quantile_vec[i]]$subscriberid],
        PlanData = CombinedMarketCounterfactual1$PlanData),
      SGpremsubs = T, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Ratio")
    risk_resbase_list[[i]][[2]]$welfare[1] <- mean(sg_surplus[alphampsi > 0 & subscriberid %in% finalsubsdata[risk_quantiles == risk_quantile_vec[i]]$subscriberid]$cs)
    risk_resbase_list[[i]][[2]]$welfare[2] <- mean(sg_surplus[alphampsi > alphampsicutoff& subscriberid %in% finalsubsdata[risk_quantiles == risk_quantile_vec[i]]$subscriberid]$cs)
    risk_res5_list[[i]] <- generateOutcomeColumns(
      CounterFactualResults = list(
        ExplodedData = CombinedMarketCounterfactual7$ExplodedData[orig_subs_id %in% finalsubsdata[risk_quantiles == risk_quantile_vec[i]]$subscriberid],
        PlanData = CombinedMarketCounterfactual7$PlanData),
      SGpremsubs = F, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Fixed")
  }
  
  ## Calculating results for family types
  ftype_resbase_list <- list()
  ftype_res5_list <- list()
  for(i in 1:4){
    print(paste("Starting family type vec", i))
    ftype_resbase_list[[i]] <- generateOutcomeColumns(
      CounterFactualResults = list(
        ExplodedData = OnlyIndividualMarketCounterfactual$ExplodedData[orig_subs_id %in% finalsubsdata[ftype == ftype_vec[i]]$subscriberid],
        PlanData = OnlyIndividualMarketCounterfactual$PlanData),
      SGpremsubs = T, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Ratio")
    ftype_resbase_list[[i]][[2]] <- generateSGOutcomeColumnsPre(
      CounterFactualResults = list(
        ExplodedData = CombinedMarketCounterfactual1$ExplodedData[orig_subs_id %in% finalsubsdata[ftype == ftype_vec[i]]$subscriberid],
        PlanData = CombinedMarketCounterfactual1$PlanData),
      SGpremsubs = T, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Ratio")
    ftype_resbase_list[[i]][[2]]$welfare[1] <- mean(sg_surplus[alphampsi > 0 & subscriberid %in% finalsubsdata[ftype == ftype_vec[i]]$subscriberid]$cs)
    ftype_resbase_list[[i]][[2]]$welfare[2] <- mean(sg_surplus[alphampsi > alphampsicutoff& subscriberid %in% finalsubsdata[ftype == ftype_vec[i]]$subscriberid]$cs)
    ftype_res5_list[[i]] <- generateOutcomeColumns(
      CounterFactualResults = list(
        ExplodedData = CombinedMarketCounterfactual7$ExplodedData[orig_subs_id %in% finalsubsdata[ftype == ftype_vec[i]]$subscriberid],
        PlanData = CombinedMarketCounterfactual7$PlanData),
      SGpremsubs = F, 
      subsidytypeInd = "Ratio", 
      subsidytypeSG = "Fixed")
  }
  
  ## Creating a Change in consumer surplus table
  foo <- data.table(
    type = c("overall", c("W. Deps, Over 50", "W.o. Deps, Over 50", "W. Deps, Under 51", "W.o. Deps, Under 51"), rep("Income Quartile", 4), rep("Risk Quartile",4)),
    value = c(1, rep(1:4, 3)),
    IM_before = rep(0, 13), 
    IM_after = rep(0, 13), 
    SG_before = rep(0, 13), 
    SG_after = rep(0, 13))
  foo$IM_before[1] <- overall_resbase_list[[1]][[1]]$welfare[2]
  foo$SG_before[1] <- overall_resbase_list[[1]][[2]]$welfare[2]
  foo$IM_after[1] <- overall_res5_list[[1]][[1]]$welfare[2]
  foo$SG_after[1] <- overall_res5_list[[1]][[2]]$welfare[2]  
  for(i in 1:4){
    foo$IM_before[1+i] <- ftype_resbase_list[[i]][[1]]$welfare[2]
    foo$SG_before[1+i] <- ftype_resbase_list[[i]][[2]]$welfare[2]
    foo$IM_after[1+i] <- ftype_res5_list[[i]][[1]]$welfare[2]
    foo$SG_after[1+i] <- ftype_res5_list[[i]][[2]]$welfare[2]  
  }
  for(i in 1:4){
    foo$IM_before[5+i] <- fpl_resbase_list[[i]][[1]]$welfare[2]
    foo$SG_before[5+i] <- fpl_resbase_list[[i]][[2]]$welfare[2]
    foo$IM_after[5+i] <- fpl_res5_list[[i]][[1]]$welfare[2]
    foo$SG_after[5+i] <- fpl_res5_list[[i]][[2]]$welfare[2]  
  }
  for(i in 1:4){
    foo$IM_before[9+i] <- risk_resbase_list[[i]][[1]]$welfare[2]
    foo$SG_before[9+i] <- risk_resbase_list[[i]][[2]]$welfare[2]
    foo$IM_after[9+i] <- risk_res5_list[[i]][[1]]$welfare[2]
    foo$SG_after[9+i] <- risk_res5_list[[i]][[2]]$welfare[2]  
  }
  foo$IM_change <- foo$IM_after - foo$IM_before
  foo$SG_change <- foo$SG_after - foo$SG_before
  bar <- foo[, .(type, value, IM_change, SG_change)]
  fred <- bar[order(type, value), ]
  
  # Saving 
  fwrite(fred, paste0(project_folder, "/analysis/tablesandfigures/release/misc/cs_change_bytype_mrkup", mrkup, ".csv"))
}